home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209b.zip / octave-2.09 / SCRIPTS.ZIP / scripts.fat / specfun / gammai.m < prev    next >
Text File  |  1997-03-08  |  3KB  |  97 lines
  1. ## Copyright (C) 1995, 1996  Kurt Hornik
  2. ## 
  3. ## This program is free software; you can redistribute it and/or modify
  4. ## it under the terms of the GNU General Public License as published by
  5. ## the Free Software Foundation; either version 2, or (at your option)
  6. ## any later version.
  7. ## 
  8. ## This program is distributed in the hope that it will be useful, but
  9. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  10. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11. ## General Public License for more details. 
  12. ## 
  13. ## You should have received a copy of the GNU General Public License
  14. ## along with this file.  If not, write to the Free Software Foundation,
  15. ## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  16.  
  17. ## usage: gammai (a, x)
  18. ##
  19. ## Computes the incomplete gamma function
  20. ##
  21. ##    gammai(a, x) 
  22. ##      = (integral from 0 to x of exp(-t) t^(a-1) dt) / gamma(a).
  23. ##
  24. ## If a is scalar, then gammai(a, x) is returned for each element of x
  25. ## and vice versa.
  26. ##
  27. ## If neither a nor x is scalar, the sizes of a and x must agree, and
  28. ## gammai is applied pointwise.
  29.  
  30. ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
  31. ## Created: 13 August 1994
  32. ## Adapted-By: jwe
  33.  
  34. function y = gammai (a, x)
  35.   
  36.   if (nargin != 2)
  37.     usage ("gammai (a, x)");
  38.   endif
  39.   
  40.   [retval, a, x] = common_s (a, x);
  41.   if (retval > 0)
  42.     error ("gammai:  a and x must be of common size or scalar");
  43.   endif
  44.   
  45.   [r, c] = size (x);
  46.   s = r * c;
  47.   x = reshape (x, 1, s);
  48.   a = reshape (a, 1, s);
  49.   y = zeros (1, s);
  50.  
  51.   k = find (!(a > 0) | isnan (x));
  52.   if any (k)
  53.     y(k) = NaN * ones (1, length (k));
  54.   endif
  55.   
  56.   k = find ((x == Inf) & (a > 0));
  57.   if any (k)
  58.     y(k) = ones (1, length (k));
  59.   endif
  60.   
  61.   ## For x < a + 1, use summation.  The below choice of L should ensure
  62.   ## that the overall error is less than eps ... 
  63.   k = find((x > 0) & (x < a + 1));
  64.   if any (k)
  65.     L = ceil (- max ([a(k), x(k)]) * log (eps));
  66.     A = cumprod ((ones (L, 1) * x(k)) ...
  67.     ./ (ones (L, 1) * a(k) + (1 : L)' * ones (1, length (k))));
  68.     y(k) = exp (-x(k) + a(k) .* log (x(k))) ...
  69.     .* (1 + sum (A)) ./ gamma (a(k) + 1);
  70.   endif
  71.  
  72.   ## For x >= a + 1, use the continued fraction.
  73.   ## Note, however, that this converges MUCH slower than the series
  74.   ## expansion for small a and x not too large!
  75.   k = find ((x >= a + 1) & (x < Inf) & (a > 0));
  76.   if any (k)
  77.     len = length (k);
  78.     u   = [zeros (1, len); ones (1, len)];
  79.     v   = [ones (1, len); x(k)];
  80.     c_old = 0;
  81.     c_new = v(1, :) ./ v(2, :);
  82.     n   = 1;
  83.     while (max (abs (c_old ./ c_new - 1)) > 10 * eps)
  84.       c_old = c_new;
  85.       u = v + u .* (ones (2, 1) * (n - a(k)));
  86.       v = u .* (ones (2, 1) * x(k)) + n * v;
  87.       c_new = v(1, :) ./ v(2, :);
  88.       n = n + 1;
  89.     endwhile
  90.     y(k) = 1 - exp (-x(k) + a(k) .* log (x(k))) .* c_new ...
  91.     ./ gamma (a(k));
  92.   endif
  93.   
  94.   y = reshape (y, r, c);
  95.  
  96. endfunction
  97.